# impulses.R
# Sims-Zha methods for computing impulse response functions (IRFs) and
# their error bands using Bayesian methods.

# Patrick Brandt
# pbrandt@utdallas.edu
# Original: August 2004
# Revised: July 2005 to use the MSBVAR library

################################################################################
# Set up functions and libraries
################################################################################
library(foreign)     # Library with functions to read Stata files.
library(MSBVAR)      # Brandt's library of BVAR functions.

################################################################################
# Read and setup data.
################################################################################
data <- read.dta("levant.weekly.79-03.dta") 
attach(data)

# Set up KEDS data as a matrix of time series.
KEDS.data <- ts(cbind(a2i,a2p,i2a,p2a,i2p,p2i),
                start=c(1979,15),
                freq=52,
                names=c("A2I","A2P","I2A","P2A","I2P","P2I"))

# Select the window or sample of data -- 1979:15--1988:50
KEDS <- window(KEDS.data, end=c(1988,50))

################################################################################
# Estimate the BVAR models 
################################################################################

# Fit a flat prior model
KEDS.BVAR.flat <- szbvar(KEDS, p=8, z=NULL, lambda0=0.6,
                    lambda1=0.1, lambda3=0.5, lambda4=0.25, lambda5=0,
                    mu5=0, mu6=0, nu=ncol(KEDS)+1, qm=4, prior=2,
                    posterior.fit=F)

KEDS.BVAR.informed <- szbvar(KEDS, p=8, z=NULL, lambda0=0.6,
                    lambda1=0.1, lambda3=0.5, lambda4=0.25, lambda5=0,
                    mu5=0, mu6=0, nu=ncol(KEDS)+1, qm=4, prior=0,
                    posterior.fit=F)

################################################################################
# Compute the impulse responses -- for flat and informed prior models
# NOTE: the mc.irf.var() mthod is generic in computing posterior
# samples of the responses.  The likelihood and Bayesian shape methods
# discussed in the Brandt and Freeman 2006 article can all be computed
# from these posterior IRF samples.
################################################################################
set.seed(116)
KEDS.irf.flat <- mc.irf.var(KEDS.BVAR.flat, 10, 5000)
KEDS.irf.informed <- mc.irf.var(KEDS.BVAR.informed, 10, 5000)

################################################################################
# Plot the impulse responses -- by various methods
# See the function documentation for the plot.mc.irf.var() function to
# see how the method argument is used to compute the difference
# responses error bands.
#
# The response plots are saved in different PDF files for viewing.
#
################################################################################

# Repsonses for the flat prior BVAR using the different error band
# construction methods.

pdf(file="IRF-percentile.pdf", width=6.25, height=8.5)
flat.impulse.percentile <- plot(KEDS.irf.flat, method=c("Percentile"), 
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Quantiles Method", outer=T, line=5)
dev.off()
        
pdf(file="IRF-normal.pdf", width=6.25, height=8.5)
flat.impulse.normal <- plot(KEDS.irf.flat, method=c("Normal Approximation"), 
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Normal Approximation Method", outer=T, line=5)
dev.off()

pdf(file="IRF-asym-eigen-1.pdf", width=6.25, height=8.5)
flat.impulse.eigen1 <- plot(KEDS.irf.flat, method=c("Sims-Zha2"), component=1,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Eigen Decomposition Quantile Method -- Component 1", outer=T, line=5)
dev.off()

pdf(file="IRF-asym-eigen-2.pdf", width=6.25, height=8.5)
flat.impulse.eigen2 <- plot(KEDS.irf.flat, method=c("Sims-Zha2"), component=2,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Eigen Decomposition Quantile Method -- Component 2", outer=T, line=5)
dev.off()

pdf(file="IRF-asym-eigen-3.pdf", width=6.25, height=8.5)
flat.impulse.eigen3 <- plot(KEDS.irf.flat, method=c("Sims-Zha2"), component=3,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Eigen Decomposition Quantile Method -- Component 3", outer=T, line=5)
dev.off()

pdf(file="IRF-sym-eigen.pdf", width=6.25, height=8.5)
flat.impulse.eigen.sym <- plot(KEDS.irf.flat, method=c("Sims-Zha1"), component=1,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Joint Normal Eigen Decomposition Method", outer=T, line=5)
dev.off()

pdf(file="IRF-asym-eigen-full.pdf", width=6.25, height=8.5)
flat.impulse.eigen.full <- plot(KEDS.irf.flat, method=c("Sims-Zha3"), component=1,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(
      main="Stacked Eigen Decomposition Method, Accounting for Correlated Repsonses",
      outer=T, line=5)
dev.off()

################################################################################
# Plot the impulse responses -- by various methods
################################################################################

# Repsonses for the reference prior BVAR using the different error band
# construction methods.

pdf(file="IRF-percentile.informed.pdf", width=6.25, height=8.5)
informed.impulse.percentile <- plot(KEDS.irf.informed, method=c("Percentile"), 
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Quantiles Method", outer=T, line=5)
dev.off()
      
pdf(file="IRF-normal.informed.pdf", width=6.25, height=8.5)
informed.impulse.normal <- plot(KEDS.irf.informed, method=c("Normal Approximation"), 
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Normal Approximation Method", outer=T, line=5)
dev.off()

pdf(file="IRF-asym-eigen-1.informed.pdf", width=6.25, height=8.5)
informed.impulse.eigen1 <- plot(KEDS.irf.informed, method=c("Sims-Zha2"), component=1,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Eigen Decomposition Quantile Method -- Component 1", outer=T, line=5)
dev.off()

pdf(file="IRF-asym-eigen-2.informed.pdf", width=6.25, height=8.5)
informed.impulse.eigen2 <- plot(KEDS.irf.informed, method=c("Sims-Zha2"), component=2,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Eigen Decomposition Quantile Method -- Component 2", outer=T, line=5)
dev.off()

pdf(file="IRF-asym-eigen-3.informed.pdf", width=6.25, height=8.5)
informed.impulse.eigen3 <- plot(KEDS.irf.informed, method=c("Sims-Zha2"), component=3,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Eigen Decomposition Quantile Method -- Component 3", outer=T, line=5)
dev.off()

pdf(file="IRF-sym-eigen.informed.pdf", width=6.25, height=8.5)
informed.impulse.eigen.sym <- plot(KEDS.irf.informed, method=c("Sims-Zha1"), component=1,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(main="Joint Normal Eigen Decomposition Method", outer=T, line=5)
dev.off()

pdf(file="IRF-asym-eigen-full.informed.pdf", width=6.25, height=8.5)
informed.impulse.eigen.full <- plot(KEDS.irf.informed, method=c("Sims-Zha3"), component=1,
            probs=c(0.025,0.16,0.84,0.975), varnames=colnames(KEDS))
title(
      main="Stacked Eigen Decomposition Method, Accounting for Correlated Repsonses",
      outer=T, line=5)
dev.off()

################################################################################
# Compute DFEVs
################################################################################
KEDS.dfev.flat <- dfev(KEDS.BVAR.flat, k=12)
KEDS.dfev.informed <- dfev(KEDS.BVAR.informed, k=12)

################################################################################
# Make graphs of just the subset of responses for I2P and P2I from
# those above.  These are the lower right 2x2 of the responses
# generated above.  Here, we combine the responses for the various
# methods in one figure for comparison.
################################################################################
# To do this, we need a specialized plot that uses the output of the
# plot.mc.irf function that can set up the series (mean, bands) for
# the impulses. Plotting needs to be handled separately. This will
# just call the various response series created above.
################################################################################

# Set up some items we need across the plots.
# Variable labels / names
varnames <- c("I2P", "P2I")

# y-axis limits
I2P.lim <- range(flat.impulse.normal$responses[,,c(29,30)],
                 flat.impulse.percentile$responses[,,c(29,30)],
                 flat.impulse.eigen.sym$responses[,,c(29,30)], 
                 flat.impulse.eigen1$responses[,,c(29,30)],
                 flat.impulse.eigen2$responses[,,c(29,30)],
                 flat.impulse.eigen3$responses[,,c(29,30)])

P2I.lim <- range(flat.impulse.normal$responses[,,c(35,36)],
                 flat.impulse.percentile$responses[,,c(35,36)],
                 flat.impulse.eigen.sym$responses[,,c(35,36)], 
                 flat.impulse.eigen1$responses[,,c(35,36)],
                 flat.impulse.eigen2$responses[,,c(35,36)],
                 flat.impulse.eigen3$responses[,,c(35,36)])

################################################################################
# Figure 1 of the Brandt and Freeman 2006 paper.
################################################################################

pdf(file="errors1-flat.pdf", width=7.5, height=5)
par(mfcol=c(2,6), mai=c(0.25,0.25,0.15,0.25), omi=c(0.15,0.75,1,0.15))
# Normal approx
ts.plot(flat.impulse.normal$responses[,c(1,4,5),29], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[1], side=2, line=3)
mtext(varnames[1], side=3, line=2)
ts.plot(flat.impulse.normal$responses[,c(1,4,5),30], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[2], side=2, line=3)
ts.plot(flat.impulse.normal$responses[,c(1,4,5),35], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[2], side=3, line=2)
ts.plot(flat.impulse.normal$responses[,c(1,4,5),36], gpars=list(xlab="",ylab=""))
abline(h=0)
# Pointwise quantiles
ts.plot(flat.impulse.percentile$responses[,c(1,4,5),29], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[1], side=3, line=2)
ts.plot(flat.impulse.percentile$responses[,c(1,4,5),30], gpars=list(xlab="",ylab=""))
abline(h=0)
ts.plot(flat.impulse.percentile$responses[,c(1,4,5),35], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[2], side=3, line=2)
ts.plot(flat.impulse.percentile$responses[,c(1,4,5),36], gpars=list(xlab="",ylab=""))
abline(h=0)
# Gaussian Linear eigenvector
ts.plot(flat.impulse.eigen.sym$responses[,c(1,4,5),29], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[1], side=3, line=2)
ts.plot(flat.impulse.eigen.sym$responses[,c(1,4,5),30], gpars=list(xlab="",ylab=""))
abline(h=0)
ts.plot(flat.impulse.eigen.sym$responses[,c(1,4,5),35], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[2], side=3, line=2)
ts.plot(flat.impulse.eigen.sym$responses[,c(1,4,5),36], gpars=list(xlab="",ylab=""))
abline(h=0)

mtext("Response in", side=2, line=3, outer=T)
mtext("Shock to", side=3, line=3, outer=T)
mtext("Normal Approximation                 Pointwise Quantiles           Normal Linear Eigenvectors", 
      side = 3, line=5, outer = T)

dev.off()

################################################################################
# Figure 2 of the Brandt and Freeman paper
################################################################################

pdf(file="errors2-flat.pdf", width=7.5, height=5)
par(mfcol=c(2,6), mai=c(0.25,0.25,0.15,0.25), omi=c(0.15,0.75,1,0.15))
# Likelihood 1
ts.plot(flat.impulse.eigen1$responses[,c(1,4,5),29], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[1], side=2, line=3)
mtext(varnames[1], side=3, line=2)
ts.plot(flat.impulse.eigen1$responses[,c(1,4,5),30], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[2], side=2, line=3)
ts.plot(flat.impulse.eigen1$responses[,c(1,4,5),35], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[2], side=3, line=2)
ts.plot(flat.impulse.eigen1$responses[,c(1,4,5),36], gpars=list(xlab="",ylab=""))
abline(h=0)
# Likelihood 2
ts.plot(flat.impulse.eigen2$responses[,c(1,4,5),29], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[1], side=3, line=2)
ts.plot(flat.impulse.eigen2$responses[,c(1,4,5),30], gpars=list(xlab="",ylab=""))
abline(h=0)
ts.plot(flat.impulse.eigen2$responses[,c(1,4,5),35], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[2], side=3, line=2)
ts.plot(flat.impulse.eigen2$responses[,c(1,4,5),36], gpars=list(xlab="",ylab=""))
abline(h=0)
# likelihood 3
ts.plot(flat.impulse.eigen3$responses[,c(1,4,5),29], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[1], side=3, line=2)
ts.plot(flat.impulse.eigen3$responses[,c(1,4,5),30], gpars=list(xlab="",ylab=""))
abline(h=0)
ts.plot(flat.impulse.eigen3$responses[,c(1,4,5),35], gpars=list(xlab="",ylab=""))
abline(h=0)
mtext(varnames[2], side=3, line=2)
ts.plot(flat.impulse.eigen3$responses[,c(1,4,5),36], gpars=list(xlab="",ylab=""))
abline(h=0)

mtext("Response in", side=2, line=3, outer=T)
mtext("Shock to", side=3, line=3, outer=T)
mtext("First Component                 Second Component                 Third Component",
      side = 3, line=5, outer = T)
dev.off()

# Print off the % of the variance explained by these 3 eigenvectors
xtable(cbind(flat.impulse.eigen1$eigenvector.fractions[c(29,30,35,36),1:5],
             rowSums(flat.impulse.eigen1$eigenvector.fractions[c(29,30,35,36), 3:5])),
       digits=rep(0,7))


# Clean up environment and save results
detach(data)
save.image(file="Levant-error-bands-analysis.RData", compress=T)

# End of file

